function SAT_SST_Step_06_post_processing(P)

    disp('Post processing starts')
    file_load = [SAT_SST_func_IO('data',P),'/Raw_SAT_anm_',P.date,'.mat'];
    load(file_load,'raw_SAT_anm')

    % Read raw temperatures ---------------------------------------------------
    file_load = [SAT_SST_func_IO('data',P),'/Raw_SAT_SST_',P.date,'.mat'];
    load(file_load,'raw_SST_anm','land_fraction','lon','lat','dis');

    % Read scaled temperatures ------------------------------------------------
    file_load  = [SAT_SST_func_IO('data',P),'/Scaled_SAT_anm_',P.date,'.mat'];
    load(file_load,'Scl_AT_anm','Scl_AT_anm_cowtan');

    % Mask SAT and SST to have the same coverage ==========================
    clear('D')
    D.raw_SAT_anm                = raw_SST_anm *0   + raw_SAT_anm;
    D.Scl_AT_anm                 = raw_SST_anm *0   + Scl_AT_anm;
    D.Scl_AT_anm_cow             = raw_SST_anm *0   + Scl_AT_anm_cowtan;
    D.raw_SST_anm                = raw_SST_anm      + raw_SAT_anm *0;
    
    % A gridded estimate of raw/scaled SAT and raw SST ====================
    disp('Calculating global and regional averages ...')
    grid_at_raw      = nan(36,18,12,171);
    grid_at_scl      = nan(36,18,12,171);
    grid_at_scl_cow  = nan(36,18,12,171);
    grid_sst_raw     = nan(36,18,12,171);
    for ct_x  = 1:36
        for ct_y = 1:18
            
            l_use = lon > (ct_x*10-10)  & lon <= (ct_x*10) & ...
                    lat > (ct_y*10-100) & lat <= (ct_y*10-90);
                
            if nnz(l_use) > 0
                % disp([ct_x ct_y])
                grid_at_raw(ct_x,ct_y,:,:)     = CDC_pnt2glbmean(D.raw_SAT_anm(l_use,:,:),   lon(l_use),lat(l_use));
                grid_at_scl(ct_x,ct_y,:,:)     = CDC_pnt2glbmean(D.Scl_AT_anm(l_use,:,:),    lon(l_use),lat(l_use));
                grid_at_scl_cow(ct_x,ct_y,:,:) = CDC_pnt2glbmean(D.Scl_AT_anm_cow(l_use,:,:),lon(l_use),lat(l_use));
                grid_sst_raw(ct_x,ct_y,:,:)    = CDC_pnt2glbmean(D.raw_SST_anm(l_use,:,:),   lon(l_use),lat(l_use));
            end
        end
    end 
    
    % Calculate trend =====================================================
    disp('Calculating trends ... ')
    clear('Trd')
    ct = 0;
    for yr_st = 1960:1:1970
        for yr_ed = 2010:1:2020
            ct   = ct + 1;
            yr   = yr_st:yr_ed;
            
            temp = CDC_trend(D.raw_SAT_anm(:,:,yr-1849),yr,3);               
            Trd.raw_SAT_anm(:,:,ct)    = temp{1};
            
            temp = CDC_trend(D.Scl_AT_anm(:,:,yr-1849),yr,3);                
            Trd.Scl_AT_anm(:,:,ct)     = temp{1};
            
            temp = CDC_trend(D.Scl_AT_anm_cow(:,:,yr-1849),yr,3);            
            Trd.Scl_AT_anm_cow(:,:,ct) = temp{1};
            
            temp = CDC_trend(D.raw_SST_anm(:,:,yr-1849),yr,3);               
            Trd.raw_SST_anm(:,:,ct)    = temp{1};
        end
    end

    % Save data ===========================================================
    file_save  = [SAT_SST_func_IO('data',P),'/Data_for_plots_',P.date,'.mat'];
    save(file_save,'grid_at_raw','grid_at_scl','grid_at_scl_cow','grid_sst_raw',...
                   'dis','land_fraction','lat','lon','Trd','D','-v7.3');
end
